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Abstract 

We present an algorithm to compute values L(s) and derivatives (s) of L- functions 
of motivic origin numerically to required accuracy. Specifically, the method applies to 
any L-series whose T-factor is of the form A s Yli=i r( s ^ 2 Aj ) with d arbitrary and complex 
Xj, not necessarily distinct. The algorithm relies on the known (or conjectural) functional 
equation for L(s). 

1 Introduction 

Many L-series in number theory and algebraic geometry can be interpreted as L-series of 
motives over number fields. For instance, Riemann and Dedekind (^-function, Dirichlet and 
Artin L-series and L-series of elliptic curves are of this kind. These are all of the form 
L(X,V,s) associated to V = H t (X) or a "motivic" subspace VcH t (X) of a projective 
algebraic variety X/K. 
Given such series 

oo 

£(s) = E^7> Res»0, (1) 

n=l 

standard conjectures state that L(s) extends to a meromorphic function on the whole of 
C and satisfies a functional equation of a predicted form. The Riemann hypothesis tells 
where the zeroes of L(s) are located and, finally, various conjectures relate values of L(s) at 
integers to arithmetic invariants of X. For instance, the Birch-Swinnerton-Dyer [||], Zagier 
|H, Deligne-Beilinson-Scholl g, || and Bloch-Kato g conjectures are examples of these. 

While the aforementioned conjectures remain unproved in the vast majority of cases, a 
lot of work has been done to provide numerical evidence for some of them in low-dimensional 



cases. This applies especially to the Riemann hypothesis for the Riemann (^-function [24|, 



Dirichlet and Artin L-series |6], 13, 15, [T^, 23] and L-series L(E, H 1 , s) of elliptic curves H 
Another well-studied case is the Birch-Swinnerton-Dyer conjecture ||, ||] for L(E,H 1 ,s)\ s =i 
where L/Q is an elliptic curve. 

To perform this kind of calculations one needs an efficient algorithm to compute L(s) 
(or, more precisely, its analytic continuation) numerically to required precision for a given 
complex s. Such algorithms are usually based on writing L(s) as a series in special functions 
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associated to the inverse Mellin transform of the T-factor of L(s). In the cases mentioned 
above these special functions are incomplete Gamma functions for dim V = 1 (Riemann 
^-function, Dirichlet characters) and incomplete Bessel functions for dim V = 2 (modular 
forms, elliptic curves). 

The original motivation for this paper was to study evidence for the "standard conjec- 
tures" in case of higher-dimensional motives, for instance curves of genus g > 1 etc. One then 
needs an efficient method to compute L(X,V,s) in cases diml/>2 as well. In this paper 
we present such a method. Namely, for an arbitrary motivic L-series for which meromorphic 
continuation and the functional equation are assumed, the algorithm numerically verifies 
the functional equation and allows to compute the values L(s) and derivatives L^ k \s) for 
arbitrary complex s to required precision. 

The scheme presented here has been implemented as a PARI script and is available on 
Q. This includes examples of computations with the Riemann C( s )> Dirichlet L-functions, 
Dedekind (^-function, Shintani's (^-function, L-series of modular forms and those associated 
to curves C/Q of genus 1,2,3 and 4. Also note that the formulae described in the paper 
can be used in any other environment as long as it provides arbitrary precision arithmetic, 
complex numbers, Laurent series and the Taylor series expansion of the T-function. 

The structure of the paper is as follows: in §|2| we start with generalities on the invariants 
of L-functions and outline the algorithm. The approach used here is standard and has been 
used in most of the algorithms to compute L-functions (e.g. [15, [19], ^]). In §|3]and §|| we 
deduce power series expansions of general Meijer G-functions required in the computations. 
In §|5] we present asymptotic expansions at infinity of the same functions and continued 
fraction expansions associated to those. Then §^ summarises the algorithm and addresses 
implementation and accuracy issues. Finally, §[?] contains some remarks on working with 
L-functions for which not all of the invariants are known. 

The author is extremely grateful to Don Zagier for suggesting to work on this project 
to begin with and also for numerous explanations, ideas and suggestions which have finally 
lead to this work. Just saying that the algorithm and this paper would not exist without his 
influence illustrates only a small fraction of the truth. The author would also like to thank 
the stimulating atmosphere of the Max-Planck Institut fur Mathematik in Bonn where most 
of this work has been carried out. 



2 Motivic L-functions 

Suppose we are given an L-series, 

oo 

L( S ) = £^, a n eC. 

n=l 

We make the following three assumptions on L(s): 

Assumption 2.1. The coefficients of L(s) grow at most polynomially in n, that is a n = 
0{n a ) for some a > 0. Equivalently the defining series for L(s) converges for Re s ^> 0. 
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Assumption 2.2. L(s) admits a meromorphic continuation to the entire complex plane. 
There exist: weight w > 0, sign e£C, real positive exponential factor A and the V -factor 

,W = r(^)...r(^) 

of dimension d > 1 and Hodge numbers Ai, . . . A^ G C, such that 

L*( S ) = A S 7 ( S )L( S ) 

satisfies a functional equatior^ 

L*(s) = eL*(w-s) . (2) 

Assumption 2.3. L*(s) has finitely many simple poles with residues rj = res s=Pj L*(s) 
and no other singularities. 

Remark. Even general motivic L-functions have much more specific parameters. Usually 
a n lie in the ring of integers of a fixed number field (most often Z), A = y/~N /ir d ^ 2 (with 
conductor NeZ), Afc are integers (or even A& G {0,1}) and e = ±l. Moreover, L*(s) is 
usually entire and there is a product formula for L(s). However, these additional assumptions 
do not simplify the algorithm and there are some L-functions not of motivic origin (e.g. 
Shintani's £- function P2|| ), to which the algorithm applies. So we do not require more than 
stated above. The last assumption that the poles of L* (s) are simple is not essential either 
(see the discussion below). 

Example 2.4. The following table contains some well-known examples of L-series satis- 
fying our assumptions and their basic invariants. For every one of those L-functions, the 
exponential factor is of the form A = ^/N lis d l 2 with N G Z. 



L(s) 


Description 


w 


d 




N 


e 


(Pi) 


C(s) 


Riemann ^-function 


1 


1 


(0) 


1 


1 


(0,1) 


L(x,s) 
L(x,s) 


X primitive Dirichlet 
character mod iV 


1 


1 


(0) , *(-!) = 1 

(1) , x(-i)=-i 


N 


|e| = 1 




C(F,s) 


Dedekind ("-function 
[F:Q]=d 


1 


d 


(o,..,o,i,...,i) 

d—a, a times 




1 


(0,1) 


Hf,s) 


f modular form 

of weight k on SL<2(Z) 


k 


2 


(0, 1) 


1 


(-l) fc 


(0,k) 


Hf,8) 


/ cusp form 

of weight k on SL2(Z) 


k 


2 


(0, 1) 


1 


(-l) fe 




Hf,8) 


/ Hecke cusp form 
of weight k on To(N) 


k 


2 


(0, i) 


N 


±1 




L(E,s) 


E/Q elliptic curve 
of conductor N 


2 


2 


(0, i) 


N 


±1 




L(C, s) 


C/Q genus g curve 
of conductor N 


2 


23 


(o,...,o,i,..,i) 

g,g times 


N 


±1 




Csh(s) 


Shintani's ("-function 


1 


4 




2 4 3 J 


1 


(o,*4,i) 



1 Functional equation may also involve two different L-functions, see Remark 2.5 
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In the second row L{x, s) satisfies a functional equation which involves the "dual" L-function 
associated to the complex conjugate character L(x, s ) ( see Remark below). In the third 
row A F is the discriminant of F j Q and a is the number of pairs of complex embeddings. 
For the last (non-motivic) example, see Shintani's original paper |22||. For the rest (and 



more), see ]18[], Chapter 4 and articles in [12] for references and additional information. For 
actual L-series computations in the above cases, see M. 



Given an L-function which satisfies |2.lH2.3] we would like to 

(a) Give a numerical verification of the functional equation for L(s). 

(b) For a given sq G C and an integer k > determine the A;-th derivative L^ k \so) to given 
precision. 

To this end define <f>(t) to be the inverse Mellin transform of 7(5), that is 

too = J o °° m t s 1 . (3) 

Henceforth we let s denote a complex number and t a positive real. The function <j>{t) exists 
(for real t>0 that is) and it decays exponentially for large t (see §|3|). In particular, the 
following sum converges exponentially fast, 

00 

e(t) = £a n <K^) (4) 



n=l 



This function is defined so that L*(s) becomes the Mellin transform of ®(t), 

n=l n=l 

roo 00 



E °» I 4>{t)O sd i = A* E J 7(-) = L*{s) . (5) 



n=l " u n=l 
By Mellin's inversion formula 

cc+ioo 



6(t) = / L*(s)t- S ds, Rec> 0, 

J c—ioo 

if c G C is chosen to lie to the right of the poles of L*(s). By the assumed functional equation 
d)forL*( S ), 

rc+ioo rc+ioo 

Q(l/t) = / L*(s)t s ds = t w eL*(w - s)t s - w ds 

J c—ioo J c—ioo 

rw—c+ioo 

= t w e L*{s)r s ds 

J w— c—ioo 
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which is almost an expression for e t w Q(t), except that the integration path lies on the left 
of the poles of L*(s). Shifting this path to the right we pick up residues of L*(s)t~ s at every 
pole of L*(s). Consequently 0(f) enjoys a functional equation, 

0(l/f) = et w e(t) - (6) 

3 

Note that the assumption that L*(s) has simple poles in not essential. If the poles are of 
higher order, the residues of L*(s)t~ s also involve some logt-terms, so (||) and ( |To|) below 
have extra terms but this does not affect the reasoning elsewhere. 

In §|| and §|sj we describe how to compute cj){t) for t > for a given T-factor 7(5). Then 
0(f) can be also effectively computed numerically since (||) converges exponentially fast. 

Now we can answer the first question, that of numerical verification of the functional 
equation for L*(s). Pick i>0 and check that (|6|) holds numerically for this t. In fact (^) 
holds for all t if and only if the functional equation (Q) is satisfied. Note that having such a 
verification is useful when not all of the invariants of L(s) are known (see §||). 

Example. Let L(s) = = Y2r^=i n ~ s be the Riemann ("-function. Then 

a n = 1, w = l, e = l, A=- T =, d = l, j(s) = T(~) . 

V 71 " 1 

We have 

00 

<t>{t) = 2e~* 2 , ®{t) = 2e-™ 2t2 . 

71=1 

The function L*(s) has simple poles in pi = and P2 = 1 with residues r\ = 1, r<i = —1, so 
the functional equation for Q(t) reads 

e(i/t) = te{t) - 1 + t. (7) 

In fact, applying Poisson's summation formula to f(x) = e~ 7Tx2 gives (0) and this proves the 
functional equation for ((s). 

We now proceed to the second problem, that of computing L(s) or L^ m \s). Fix s G C 
and let 

G s (t) =t~ s I™ <j){x)x s — , t>0. (8) 
Jt x 

Thus t s G s (t) is the incomplete Mellin transform of <p(t) and lim^o t s G s {t) = 7(5) is the 
original T— factor. Again the function G s (t) decays exponentially with t and can be effectively 
computed numerically (§|4],||). 

Consider (||) which expresses L*(s) as the Mellin transform of 0(t). Split the integral 
into two and apply functional equation (|6|) to the second integral: 

l*( S ) = r e(t)t s - = r + f 1 = r 9(0^- + r e(i/t)t- s - = 

Jo 1 Ji Jo Jl t Jl 1 

^ @(t)t s ^ + J™ et w Q{t)r sd ± - /°° E r^r sd l = (9) 

/°°e(t)t s - + e re(t)t w - s ^ + Y^^. 

Jl t J x * jPj-s 
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By definition of @(t) and that of G s (x), the first integral can be rewritten, 

e(t)ff= / E^(f)* s f = £^/ ^ s f = 

n=l n=l 
n=l n * n= i A 

The same applies to the second integral (replace s by w — s), so @ becomes 

oo oo 

L*( S ) = Y: a n G s q) + e £ a n G._ s (5) + E ^ • 

n=l n=l j 

This formula allows to determine L*(s) and hence L(s)=L*(s)/'y(s) for a given s G C. 
Differentiation also gives the formula for the derivatives, 

, oo , oc , 

h L 'M = E -& G -<7> + E E ^-.(3) + E ff^ do) 

n=l n=l j 

It remains to explain how to compute the functions <f>(t) and -j^G s (t). This occupies the 
next three sections. 

Remark 2.5. We assumed that the functional equation (|2|) involves L*(s) both on the left- 
hand and on the right-hand side. In fact, for general motives the functional equation may 
be of the form 

L*(s) = eL*(w-s) 

where 

00 00 ^ 

L ( a ) = y^, L(s) = y?n 
w ^ n s w ^ n s 

?1=1 71=1 

are L-functions of "dual" motives. For instance, Dirichlet L-series associated to non-quadratic 
characters are of this nature. Clearly our arguments go through in this more general case as 
well. The result is that @ and fllQ| ) have to be simply replaced by 

e(i/t) = et w e(t) - ^f/i. 

and 

S^-m = E -^a.® + £ E + E 

n=l ra=l j 

Here -A,j5j- etc. are associated to L(s) as etc. are to L(s). 
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3 Computing <f)(t) for t small 

Recall that 

7W = r( £^)...r(£^) ,n, 

and (ft(t) is denned as the inverse Mellin transform of 7(5). By Merlin's inversion formula 
(see e.g. J3J, §2), (f)(t) is given by a residue sum, 

4>{t) = J2 res ^ l(.s)r s , i >0 . (12) 
zee 

Since T(s) has simple poles at and the negative integers, the function 7(5) has a pole at 
s 6 C iff s = —Xj — 2n for some j and an integer n. If Xj — A& 2Z for j 7^ k, then all poles 
of 7(s) are simple and 



res s 



- Aj - 2 „(7(sr s ) = 2 { -^ft^ n.7( ( 



Mi 



Hence in this case (12) is of the form J^j^Pji^ 2 ) where Pj(t) a power series in t. The 
coefficients of pj (t) satisfy a simple linear recursion coming from the relation T(s + 1) = sT(s) . 

Example 3.1. Let d = 1 and let Ai be arbitrary. Then <p(t) is given by 

00 ( —-\\ n 
<j>(t) = t A i y 2 ^^t 2n = 2i Al e-* . 

t'o n! 

In general, the poles of 7(s) are not simple and the residue of ^y(s)t~ s in s = z is t~ z 
times a polynomial in hit of a fixed degree. The reason is that non-constant terms of the 
Taylor expansion of t~ s around s = z, 



Int 



fc=o K - 



(s — z 



contribute to the residue in case of a multiple pole. So ( |12| ) is again of the form J2j t ■ Pj(t • )> 
except now is a power series in i whose coefficients are polynomials in lni of a fixed 
degree depending on j. 

Example 3.2. Let d = 2 and Ai = A2 = 0. Then (p(t) is a Bessel function, 

<f>{t) = 4K (2t) = -4(lnt+ 7e ) - 4(lni-l+ 7e )t 2 - 21nf ' 3+2 > t 4 + . . . 
where 7 e = — r'(l) is the Euler constant. 

Algorithm 3.3. Expansion of <p(t) for t small. The following describes the recursions 
necessary to determine the coefficients of (^) for a general T-factor 7(5). 
1. Let 7(s) and 4>(t) be defined by ( |i~l| ) and ( |i"2j ) respectively. 
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2. We say that Xj and A& are equivalent if Xj — Afc£2Z. Let Ai,...,Ajv denote the 
equivalence classes and let lj = \Aj \ . Thus J2 h = d- 

3. Let rrij = — A&. + 2 where A& . G Aj is the element with the smallest real part, that is 
inf^gAj Re A = Re A^ . In other words, 7(5) is analytic at s = rrij, has a pole (of some order) 
at s = rrij — 2 and a pole of order lj at s = rrij — 2n for n 3> 0. 

4. Let cj 0) (s) be the be ginning of the Taylor series of 7(5 + mj) around s = which ends 
with 0(s lj ) as the last term. 

5. For 1 < j < d and n > 1 define Cj(s) recursively by 

c f\s) = cf- 1 \s)/l[( s -±^-n), (13) 

fc=i 

considered quotient of Laurent series in s = 0. Note that cf\s) ends with 0(1) for 

n»fl. Let c^l denote the coefficient of s~ k in c^ n ^(s). 

6. For t real positive, (f>(t) is given by 

N 00 Jj-1 

#) = E^E(E^&)* 2n - ( 14 ) 

j=l n =l fc=0 

Remark 3.4. The series above converges exponentially fast since 

max Ic^J = 0((n\)~ d ), as n — > 00 . 

j<N,k<lj h ~ k 

However, for large t this is not a very effective way to compute (f>(t). Take for instance the 
series e~* = X)n^=o( — t) n /n\ for i = 50. The terms grow up to 3. x 10 20 for n = 50 before 
starting to decrease to in absolute value. Thus to determine e -50 to 10 decimal digits with 
this series one has to require working precision of 30 digits and compute 160 terms until 
everything happily cancels leading to the answer 0.0000000000. This is clearly not terribly 
effective. As this is exactly the general behaviour for arbitrary 7(5) we use a different method 
for large t, based on the asymptotic expansions at infinity. This is described in §[| below. 



4 Computing G s (t) for t small 

As explained in §^, we also need a way to compute the incomplete Mellin transform of 4>(t) 
and its derivatives. Recall that for s € C and t > we define G s (t) to be 

roc a t 

G s (t)=t~ s / <p(x)x s — . 

Jt x 

Recall also that limj^o t~ s G s (t) exists and equals 7(5) whenever s is not a pole of j(s). For 
such s clearly 

t s G s (t) = 7 ( S ) - f^x 8 — . (15) 
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Since (|TJ) expresses 4>(t) as an infinite sum of terms of the form t a (lnt)P, term-wise integra- 
tion in (p3) gives a similar expression for G s (t). 

In the points where 7(5) does have a pole, the formula ( |l5|) does not make sense as the 
right-hand side becomes 00 — 00. However, it is not difficult to locate the terms which give 
contributions to the principal parts of the Laurent series. Then ignoring these terms gives 
the value of G s (t) for such s. 



Algorithm 4.1. Expansion of j^Gsft) for t small. This is all summarised in the 

Qk 



following formulae which allow to determine -^-^G s (t) for arbitrary s £ C and t > 0. Here 



a E C and i,j,k>0,n>l are integers. 
1. Let c { r) be as in (111). 



2. Define L a j^(x) & C[x] by the formula 



L a , hk {x) = { -^=0 V * 



a = 0. 



3. Let 



:( n ) ^ _ „(«) 
i=l 

4. For i > consider the infinite sum 

N 



S tls( X ) = £ 4" ) L 2 n+ S -m ;7 ,i,fc(x) G C[z] 



G.,*(t) = £ ^ E S tU^) t 2n (16) 
i=i «=i 

5. The formula for 4^G s (t) is 

where f(S)* s=s denotes the constant term clq of the Laurent expansion J2k a k{S — s) k of f(S) 
at S = s. Thus f (S)* s=s = f (s) if f(S) is analytic at S = s. 



Remark 4.2. The series for j^G s (t) converges exponentially fast since the corresponding 
one for <p(t) does (cf. 3.4). Again, however, it is not effective for large t, in which case we 
use an alternative approach described in the following section. 



5 Computing <p(t) and G s {t) for t large 

To compute 4>(t) and G s (t) for large t we begin with the asymptotic expansions for these 
functions around infinity. 

Recall that <p(t) is defined as the inverse Mellin transform of a product of T-functions, 

r(^)--(^H>K?. 



9 



In other words, (f)(t) is a special case of Meijer G-function. Given two sequences of complex 
parameters, 

oil • • • > o, n , a n +i, . . . a p and bx, b m , & m +i> ... b q 
a general Meijer G-function G™^ n (i; ax, a p ; bi, b q ) is denned by 



oo 

s 



dt \-[™ =l T{b j+ s)Y$ =l T{l- 



We refer to Luke jL7]], 5.2-5.11 for basic properties of the G-function. 
In our case replacing s by s/2 yields an identification 

m = 2 Gg(t 2 ; ; ^) 

As discovered by Meijer (in greater generality), the function Gq'^ possesses the following 
asymptotic expansion at infinity ( |i~7|| , Theorem 5.7.5) 

oo 
ii.=0 

d 

« = (l-d+^A i )/2 (17) 
i=i 

Here M n = M n (Ai, A^) are constants, Mo = 1. As for 4>(t), it follows that 

oo 

0(^/2) ^ 2(2 ^ e-^t^Mnr" (18) 

n=0 

We would like to note here that the stated asymptotic expansion for large t is much 
"neater" than the expansion ( |l4|) of 0(t) for small t: it does not involve any logarithmic 
terms and its shape is independent of whether any of the \j are equal modulo 2Z. 

The coefficients M n in the asymptotic expansion can be found as follows. The defining 
relation 



7 (s + 2) = 7 («) x J ; 



s + Xj 



on the level of inverse Mellin transforms is equivalent to an ordinary differential equation 
(of degree d) with polynomial coefficients for 4>{t). It follows that the function t~ K e dt (j)(t d / 2 ) 
satisfies a different ODE, of degree d+1. Formally substituting 1 + J2 n >i M n t~ n as a 
solution gives a recursion for the M n with polynomial coefficients. This has been worked 
out in general by E. M. Wright; see Luke |i~7| ], 5.11.5, especially formulae (8) and (16) for 
details. 

Algorithm 5.1. Asymptotic expansion associated to 4>(t). Here is the answer in our 
case, re-written in a slightly different polynomial basis. 
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1. Let S m = S m (Xi, Ad) denote the m-th elementary function of Ai, Ad, 



d d 

So = 1, Si = E A 7 -, . . . , Sd = JJ Aj . 

j=i i=i 

2. Define also modified symmetric functions Sd+i = and 

m 

S m = , Z{-Si) k d m - 1 - k ( k+ i- m )S m -h, 0<m<d. 

k=0 

3. For k > construct Ak(x) G Q[x] by means of the generating function 



,2k 



k=0 



4. For p > 1 consider the following polynomials 



u p (n) 



d 



1 i p~" 

? E^IK d -i) E 



(2n-p+l)P- m - 2fc 



('2d)?' ^-^ J-- 1 - v " " ' ' ^— ' (p— m— 2fc)! 
^ ' m=0 i=»ri fc=0 



A fc (d-p) 



5. The coefficients M n in the asymptotic expansion (|18|) satisfy a recursion 
M n = < 



0, n < 0, 

1, 71 = 0, 

1 w 



I ~ Ep=i ^(n)M n _ p , n > 1. 

Applying term- wise integration to (fiq), it is also easy to deduce the asymptotic expansion 
for G s (t) for t— >oo, 



G s (i d / 2 ) 



(19) 



n=0 



Here /c = (1 — cZ + Si)/ 2 as in (O) and /^ n (s) = /^ n (Ai, Ad! s) satisfy a recursion 



0, 
1, 

p=l 



n 



Si+rf(s-l)-2(n-p)-l 
2d 



n < 0, 
n = 0, 

v p {n))n n - p , n>l. 



(20) 



By induction one shows that /% is a polynomial in s with the leading term 2 n s n . So if we 
differentiate ( |I9| ) A; times to s, exactly terms vanish and we get the following formula for 



the derivatives J^G s (i), 



ds k 



G s {t d > 2 ) 



(2 7 r)( d - 1 >/ 2 e -di t K-l-k 9 fc Mn+fc(^) r 



n=0 



(21) 
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Equations ([lS|), ( |l9| ) and (21) provide asymptotic series for the functions (j>(t), G s (t) and 



JprG s (t) at infinity. For computational purposes, though, it is better to work with continued 
fraction expansions associated to these series. Consider, for instance, the case of 4>{t), the 
case of J^G s (i) being analogous. 

Fix d and Ai, A^. Letting x = 1/t in (|l8|) we get 

_ 00 



ra=0 



with M n constants. As any formal series, the right-hand side can be formally written either 
as a unique infinite continued fraction 



J2 M nX n = a + — , a n ^0 for n>0, 

n=o «i + 



*3 + -- 



or as a unique terminating fraction of the same form. To see this start with po(x) = M n x n 
and define recursively formal power series p n +\(x) in terms of p n (%) by 

p n {x) =a n -\ —-, n>0. 

Pn+l(X) 

Here k n is the degree of the first non-zero term in p n (x) — Pn(0); if Pn( x ) — f° r some n, 
then terminate. This shows the existence of the continued fraction expansion and uniqueness 
is not difficult to verify as well. This construction also gives an explicit way to calculate the 
a n if the M n are given (up to some limit). There are of course better (computationally more 
stable) methods, see for instance fll|, [Tgl - 

If the fraction does not terminate, define the partial convergents C n {x) for all n by 



C n {x) = a + 



«1 + 



If the fraction does terminate at Cat, let C n = Cm for n> N . 

In any case we can think of C n {x) as approximants to the original function \{j(x) of (p^). 
Indeed, C n {x) is a rational function whose Taylor expansion around x = starts at least 
with Mq + ... + M n x n . Hence il)(x) and C n {x) have the same asymptotic expansions up to 
x n . Thus there is a constant K n > such that 

\*P(x) - C n (x)\ <K n x n+1 . 

Unfortunately, it seems very difficult to provide explicit bounds for K n . It appears that 
C n (x) converge rapidly to ip(x) but to prove either "converge" or "rapidly" or "to ip(x)" in 
any generality seems hard. So the last step of the algorithm is based purely on empirical 
observations concerning the convergence of the continued fractions. If one is uncomfortable 
with this, see §^ on how to avoid this. In the implementation Q we do use asymptotic 
expansions with a simple numerical check (see step 7 below) to justify the values. 
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Algorithm 5.2. Computing 0(t) for f arbitrary. The computation of 0(f) for arbitrary 
t can be performed as follows: 

1. Let e > be the necessary upper bound for the required precision in the computations 

of <t>(t). 

2. Let n (f) be the n-th approximant to 0(t) defined by (cf. (|l8|)) 

Mt) = 2(2,)^ e - dt ^ t 2./ d c^l/t^), n > . 

As we already noted, 0(f) — 4>n(t) = 0(t~ n ) as f — » oo . Denote by (j>tayior(t) the function 0(f) 
computed using the power series expansion at the origin as in §[|. 

3. Determine to such that |0o(t)| < e/2 for t > to- 

4. Choose a subdivision of the interval [0, to]) 

< tfe < tk-i < • • • < t\ < to < oo 

For every tj let rij be an integer for which \(j> n (t) — n+ i(t)| < e/2 and \4> n (t) — n +2(i)| < e/2. 

5. Determine M n for < n < rife using the recursion they satisfy. 

6. The function 0(f) is computed as follows 

4>taylor(t), < t < tfc 
4>general(t) = ^ 0n;(f), tj < t < fj-1 

0, t > t 

7. As a numerical check, verify that 4>taylor(tk) agrees with 0n fe (tfc)- 



Example 5.3. Let d = 2 and Ai = A2 = as in Example 3J2. Recall that 0(f) = 4Ko(2t) 
is a Bessel function in this case. Asymptotic expansion (|l8[) then reads 



0(f) ~ 2v^F e- 2 * t~ 1/2 ^ M n t" n 

n=0 

and the coefficients M n satisfy a recursion 

16nM„ = -(2n - l) 2 M n _i . 

It follows that 

M = l, M 1 = -i M 2 = -g^, ■■■,M w = - ^- 1 B- 1) " ,... 

Choose e=5 ■ 10~ 10 and to = 12, fx = 6, ti = 2. Take n\ = 6 and n2 = 20 and compute 
^(f) by 

<t>taylor(t), < t < 2 



620 (t), 2 < t < 7 

6 (f), 7 < t < 12 

0, t > 12 



4* general^f) - 

[ 0, t > 12 

As a numerical check we verify that \4>tayior(2) — 02o(2)| < 4 • 10 -14 < e as required. 
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6 Implementation notes 



Let us begin with a summary of steps of the algorithm presented in the previous sections. 
We start with an L-function satisfying [2,1| - |2,3| (see also Remark p,5| ). 



The formula used for the numerical verification of the functional equation is @ and 
that for computing L(s) and its derivatives is ( |ToD together with L*(s) = L(s)/ , y(s). 
The functions used in these formulae are <fi(t) defined by @ and G s (t) defined by @. 



To compute (j)(t) numerically we employ ( [5,21 ). It is based on power series expansions in 
the origin (|3.3|) , asymptotic formula (|l8|), recursion ( |5.1|) and the associated continued 
fractions. 

Similarly J^G s (i) is computed in the same way. The corresponding expansion in the 



origin is given by (4.1), asymptotics by (W% and recursion for the coefficients by (pQ). 



However, in order to make a practical algorithm out of these results, we still need to 
explain how to truncate various infinite sums and discuss related precision issues. 

If one desires to implement our method with rigorous proofs that all of the computations 
are correct, the following issues have to be considered. First, one has to keep track of the 
number of operations used and the possible round-off errors, perhaps even using interval 
arithmetic to justify the computations. Second, one needs to have analytic bounds on the 
size of the functions <j){t) and G s (t) for large t, rather than just asymptotic behaviour. 

In the PARI implementation Q we have chosen to be content with the intuitively natural 
bounds and a few numerical checks to justify the results. A reader wishing to use a more 
rigorous approach, might consider the following: 

Remark 6.1. Let us start with the computations of <p(t) and -^G s (t) by means of series 
expansions around the origin. Both ( |l4] ) and (|i~6| ) defining these are infinite sums, but it is 
not difficult to see how to terminate them. The point is that it suffices to give an explicit 



bound on the coefficients cf^ which goes to exponentially with n. Everything else in (14) 



and (|lj) grows at most polynomially in n so any rough estimate on them will do. As for an 
explicit exponential bound on , it follows from (|l3|) and, say, the obvious lower bound 
^n d for n>0 on the coefficients in ITfc=i (1 — S+ 'i +mi ) treated as a polynomial in s. 
Remark 6.2. The next question is that of working precision. This has already been men- 



tioned in 3.4 and |4.2| . When <fi(t) and G s (t) are computed from the expansions around the 
origin, the terms in the defining series can be very large. Therefore one needs to work with 
larger working precision than the desired precision of the answer. 

Similar cancellation in fact occurs when one computes L(s) for s with large imaginary 
part. This is a well-known problem, which one has to face when verifying the Riemann 
hypothesis. The point is that L(s) = L*(s)/'j(s) and both L*(s) and 7(s) decrease expo- 
nentially fast on the vertical strips. Hence one needs to compute L*(s) to more significant 
digits (log 10 |7(s)| more to be precise), to evaluate L(s) to given precision. 

A solution to this has been suggested in Lagarias-Odlyzko |l5| and worked out by Rubin- 



stein ]19[, at least for d = 1 and d = 2. By modifying G s (t) by a suitably chosen exponential 
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factor one obtains a formula for L(s) which does not have the loss-of-precision behaviour. 
It is perhaps possible to use Rubinstein's approach to obtain similar formulae for a general 
T— factor as well. 

Remark 6.3. The next issue is that of truncating the main formulae used is this paper. 
Recall that to verify the functional equation numerically, we use the function 

oo 

e(*) = X>n<Kf)- (23) 

71=1 

Then to actually compute the L- values, we use 

= E ^ k G a Q + e £ a n ±- k G^ s q) + £ ^§ . (24) 

n=l n=l j 



See also Remark 2.5 for the necessary modifications when there are two different L- functions 
involved in the functional equation. In any case, one needs analytic estimates on <f>(t) and 
-7^G s (t) for large t to carefully estimate the error in truncating these infinite sums. 

One possible way to obtain such estimates is to use Tollis' method |23| based on Braaksma's 
work |H on asymptotic behaviour of Meijer G-functions. By applying the Euler-Maclaurin 
summation formula to the Merlin-Barnes integral defining G s (t), Tollis determines an explicit 
exponential bound for G s (t) in the case A = (0, 0, 1, 1) with p + a zeroes and a ones 
(p, cr>0). It is likely that his method is general enough to obtain similar estimates for an 
arbitrary T-factor as well. 

Remark 6.4. Finally, let us turn to the methods of §||, asymptotic expansions and associ- 
ated continued fractions. 

Unfortunately, there seem to be few cases where one can actually provide explicit esti- 
mates for the convergence of the continued fractions of, say, <p(t). The most general result 
known to the author in this connection is that of Gargantini and Henrici (TO]. They show 
that functions which can be written as Stieltjes transforms of positive measures admit con- 
vergent continued fraction expansions at infinity, with explicit error bounds. This does 
not seem to apply to our functions in general, though. See Henrici |11], Chapter 12 and 



Lorentzen-Waadeland 16] for more information. 

The analysis is available, though, in low-dimensional cases. For instance for d = 1 the 
function G s (x) is the incomplete Gamma function, for which there are known convergent 



continued fractions expansions at infinity, see Henrici [11], 12.13.1. Also for d = 2 the function 
4>{x) reduces to the Whittaker function, which is a Stieltjes transform (basically, of itself). 
So in this case the continued fraction expansion converges, see Henrici Q, 12.13.11. 

One possible way out is to compute 4>{t) and -^G s (t) only using Taylor expansions at the 
origin, even for large t. It is easier to give precise estimates for the convergence in this case, 
although one does pay the price with substantial loss of efficiency Alternatively, one might 
try a completely different approach to compute the functions in question at infinity. For 
instance, it is perhaps possible to use backward recursions, as one does for Bessel functions. 
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7 L-functions with unknown invariants 



In a perfect world, one knows all of the invariants associated to one's L-function. In a less 
perfect world, one does not know exactly the sign e and, perhaps, the residues at the poles 
of L*(s). In reality, however, there are plenty of examples where it is difficult to determine 
the exponential factor A and even some of the coefficients a n . Fortunately, in some of these 
cases it is still possible to make computations with L-functions. 



To illustrate this, say that L(s) is expected to satisfy Assumptions 2.1 - |2.2| and only the 
sign e in the functional equation is difficult to determine. As we already mentioned, the 
functional equation is equivalent to the statement that for all 1 < t < oo , 

6(1/*) = et w O(t) - Y, r i tPj ■ ( 25 ) 

3 

Choose 1 < t < oo and evaluate the left-hand and the right-hand side. This gives an equation 
which can be solved for e. Of course it is then sensible to verify that fl25| ) holds with the 
obtained e by verifying it numerically for some other values of t. 

The same applies in the case when neither e nor the residues ri are known. The equation 
above is linear in all of these, so choosing enough t's gives a linear system of equations from 
which e and the rj can be obtained. There might be of course precision problems if there are 
many residues to be determined. 

In most cases, actually, e = ±1 and L*(s) has no poles, so simply trying e = —1 and e = 1 
for some t > 1 immediately yields the right sign. 

Next come the dimension d, the Hodge numbers Ai,...,A<2 and the poles pj of L*(s). 
Fortunately, these can always be determined in practice, at least in all of the cases that the 
author is aware of. 

The next issue is that of the exponential factor A. For instance, consider L(C, H 1 , s), 
the L-function associated to H 1 of a genus g curve C/Q. Then A = y/N '/ir 9 where N 
is the conductor of C. In practice, to determine N one needs at least to be able to find 
a model of C over Z which is regular at a given prime of bad reduction. This, in turn, 
means performing successive blowing-ups over the unramified closure of Q p , an operation not 
without computational difficulties. For curves of genus 1 and 2 there are effective algorithms 
for doing this, but not for higher genus. So finding N for a given curve might be hard in 
practice. Also note that (25) is absolutely not linear in A, so one cannot solve for it directly. 



Fortunately, one can usually determine the full set S = {pi, •• of primes where C 
has bad reduction. Then one knows that N = p bl ■ ■ -p b k k is composed of those primes and 
has (hopefully) an upper bound for the 6j, say in terms of the discriminant of C or some 
similar quantity. This leaves only finitely many choices for N and (as in the case of the sign 
e = ±1), a simple trial-and-error can establish the proper functional equation. It should be 
noted here that this applies, of course, only to those L-functions for which there is a unique 
A (and e etc.) for which the functional equation holds. 

Finally we come to the coefficients a^. Again take the case of a genus g curve C/Q with 
the set E = {pi, ...,pk} of bad primes as above. Then the problem is to determine the local 
factors at bad primes, that is the coefficients a pJ for 1 < i < k and j > 1 . The local factors 
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at good primes can be determined by counting points over finite fields and the coefficients 
a n for composite n can be obtained by the product formula. 

Fortunately, again, there are only finitely many choices of possible local factors for a 
given bad prime pi. For instance \a p \ < 2 A 7p and \a p j\ satisfy similar estimates. Moreover, 
a p j with 1 < j < 2g determine a p] for all j as the degree of the local factor is bounded by 
2g. Note, however, that this is not a very practical approach, especially for large primes pi 
when there are numerous possibilities for the local factors. 

Another approach is to note that the functional equation (^) is in fact linear in the dj, 
since Q(t) is. If there were only finitely many unknown coefficients aj, they could obtained 
in the same way as e and the rj were. 

To illustrate what can be done when infinitely many coefficients are unknown, consider 
the following typical case: 

1. Say, there is just one prime p for which a p is difficult to determine theoretically, 

2. assume that all a n are integers, 

3. assume that there is a product formula for L(s) in question, in particular a mn = a m a n 
for m, n co-prime. 

Using multiplicativity, write @ as 

oo 
k=l 

where 9k(t) are computable functions. Moreover, since we only take finitely many terms 
when actually computing Q(t), we have 

K 

©(*) ~ E a p^k(t) 

k=l 

where "~" stands for "equal to required precision" . Hence the functional equation ( p5| ) for a 
fixed t becomes simply a linear equation in a p , a p K. Thus we can again plug in enough t's 
to get a linear system which can be solved for the a p %. However, the coefficient functions 6k(t) 
decay rapidly with k, so a p % obtained from solving this system are certainly unreliable for 
large i. If, however, the first coefficient a p does look like an integer, we can simply round it 
off and repeat the same process with a p 2, a p n as variables until all the a p % are determined. 

In practice this works well for a large prime p and even when there are several (large) 
primes p for which the ay are unknown. This does not work for small primes, for instance 
virtually never for p = 2. But then for small primes one may try all possible local factors 
by trial-end-error and for large primes solve for the coefficients as described here. 

At this moment the reader might be long horrified by the methods suggested here and 
might wonder whether the reliability of such an approach is not extremely dubious. In our 
defence we can say that since there is a very effective way to verify the functional equation 
numerically, any method to make an intelligent guess will do, however dubious it might be. 
When A, e and the bad local factors are determined (or simply guessed in whatever way), one 
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can make numerous checks that these have the correct shape and that (|6|) holds for various 
t. Thus it is hoped that someone who actually tried to perform blowing-ups on a genus 
6 arithmetic surface which has 2 20 in the discriminant will forgive the author for offering 
desperate tricks to avoid the hard work. After all, this does allow to give evidence for various 
conjectures even in the difficult cases where it is hard to determine all of the invariants of 
the L-function in question using theoretical arguments. 

Finally, let us mention here that there are fortunately better ways to guess the local 
factors for bad primes, at least for arithmetic surfaces. These have been used to make 
computations with curves of genus g < 8 and are to appear in ||. 
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